clear all

set more off
set matsize 10000

global base "X:\Dropbox\educacion superior 2021\submission\EDCC\For publication\Replication package\Data"
global output "X:\Dropbox\educacion superior 2021\submission\EDCC\For publication\Replication package\Output"

cd "$base"
use "output/panel_final"

*===============================================================================
* Main regressions
*===============================================================================

*** 1. Only wealth (T3 - A)

eststo main_1: sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(terwi2_r3 terwi3_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model
{
	
cap mat drop ModelA_main
mat ModelA_main = J(26,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 {
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat ModelA_main[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat ModelA_main[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}

}

** exploring the wealth gradient
{
	
cap mat drop ModelA_wealth
mat ModelA_wealth = J(10,2,.)

forvalues s=0/1 { // male then female
	
	* bottom vs middle tercile
	
	qui lincom _b[enroll_sup:`s'.sex#c.terwi2_r3]
	local b_`s' = `=round(`r(estimate)', 0.01)'
	local se_`s' = `=round(`r(se)', 0.01)'
	local pval_`s' = `r(p)'
	
	local c = 1
	foreach coef in b se { // coef then se
		mat ModelA_wealth[`=`s'*2+`c'',1] = ``coef'_`s''
		local c = 2
	}

	mat ModelA_wealth[`=`s'*2+1',2] = `=cond(`pval_`s''<0.01, 3, `=cond(`pval_`s''<0.05, 2, `=cond(`pval_`s''<0.1, 1, 0)')')'

	* middle vs top tercile

	qui lincom (_b[enroll_sup:`s'.sex#c.terwi3_r3] - _b[enroll_sup:`s'.sex#c.terwi2_r3])
	local b_`s' = `=round(`r(estimate)', 0.01)'
	local se_`s' = `=round(`r(se)', 0.01)'
	local pval_`s' = `r(p)'
	
	local c = 1
	foreach coef in b se { // coef then se
		mat ModelA_wealth[`=5+`s'*2+`c'',1] = ``coef'_`s''
		local c = 2
	}

	mat ModelA_wealth[`=5+`s'*2+1',2] = `=cond(`pval_`s''<0.01, 3, `=cond(`pval_`s''<0.05, 2, `=cond(`pval_`s''<0.1, 1, 0)')')'

}

qui test _b[enroll_sup:1.sex#c.terwi2_r3] = _b[enroll_sup:0.sex#c.terwi2_r3]
mat ModelA_wealth[5,1] = `=round(`r(p)', 0.001)'

qui test (_b[enroll_sup:1.sex#c.terwi3_r3] - _b[enroll_sup:1.sex#c.terwi2_r3]) = (_b[enroll_sup:0.sex#c.terwi3_r3] - _b[enroll_sup:0.sex#c.terwi2_r3])
mat ModelA_wealth[10,1] = `=round(`r(p)', 0.001)'

}

*** 2. Family structure (T3 - B)

eststo main_2: sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 migrated_r1 indigena_r1 order_1 msib_r3 fsib_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model
{
cap mat drop ModelB_main
mat ModelB_main = J(26,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 migrated_r1 indigena_r1 order_1 msib_r3 fsib_r3  {
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat ModelB_main[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat ModelB_main[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

** exploring the wealth gradient
{
	
cap mat drop ModelB_wealth
mat ModelB_wealth = J(10,2,.)

forvalues s=0/1 { // male then female
	
	* bottom vs middle tercile
	
	qui lincom _b[enroll_sup:`s'.sex#c.terwi2_r3]
	local b_`s' = `=round(`r(estimate)', 0.01)'
	local se_`s' = `=round(`r(se)', 0.01)'
	local pval_`s' = `r(p)'
	
	local c = 1
	foreach coef in b se { // coef then se
		mat ModelB_wealth[`=`s'*2+`c'',1] = ``coef'_`s''
		local c = 2
	}

	mat ModelB_wealth[`=`s'*2+1',2] = `=cond(`pval_`s''<0.01, 3, `=cond(`pval_`s''<0.05, 2, `=cond(`pval_`s''<0.1, 1, 0)')')'

	* middle vs top tercile

	qui lincom (_b[enroll_sup:`s'.sex#c.terwi3_r3] - _b[enroll_sup:`s'.sex#c.terwi2_r3])
	local b_`s' = `=round(`r(estimate)', 0.01)'
	local se_`s' = `=round(`r(se)', 0.01)'
	local pval_`s' = `r(p)'
	
	local c = 1
	foreach coef in b se { // coef then se
		mat ModelB_wealth[`=5+`s'*2+`c'',1] = ``coef'_`s''
		local c = 2
	}

	mat ModelB_wealth[`=5+`s'*2+1',2] = `=cond(`pval_`s''<0.01, 3, `=cond(`pval_`s''<0.05, 2, `=cond(`pval_`s''<0.1, 1, 0)')')'

}

qui test _b[enroll_sup:1.sex#c.terwi2_r3] = _b[enroll_sup:0.sex#c.terwi2_r3]
mat ModelB_wealth[5,1] = `=round(`r(p)', 0.001)'

qui test (_b[enroll_sup:1.sex#c.terwi3_r3] - _b[enroll_sup:1.sex#c.terwi2_r3]) = (_b[enroll_sup:0.sex#c.terwi3_r3] - _b[enroll_sup:0.sex#c.terwi2_r3])
mat ModelB_wealth[10,1] = `=round(`r(p)', 0.001)'

}

*** 3. Skills (T3 - C)

eststo main_3: sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 migrated_r1 indigena_r1 order_1 msib_r3 fsib_r3 skills_cr_r3 skills_nc2_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model
{
cap mat drop ModelC_main
mat ModelC_main = J(26,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 migrated_r1 indigena_r1 order_1 msib_r3 fsib_r3 skills_cr_r3 skills_nc2_r3 {
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat ModelC_main[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat ModelC_main[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

** exploring the wealth gradient
{
	
cap mat drop ModelC_wealth
mat ModelC_wealth = J(10,2,.)

forvalues s=0/1 { // male then female
	
	* bottom vs middle tercile
	
	qui lincom _b[enroll_sup:`s'.sex#c.terwi2_r3]
	local b_`s' = `=round(`r(estimate)', 0.01)'
	local se_`s' = `=round(`r(se)', 0.01)'
	local pval_`s' = `r(p)'
	
	local c = 1
	foreach coef in b se { // coef then se
		mat ModelC_wealth[`=`s'*2+`c'',1] = ``coef'_`s''
		local c = 2
	}

	mat ModelC_wealth[`=`s'*2+1',2] = `=cond(`pval_`s''<0.01, 3, `=cond(`pval_`s''<0.05, 2, `=cond(`pval_`s''<0.1, 1, 0)')')'

	* middle vs top tercile

	qui lincom (_b[enroll_sup:`s'.sex#c.terwi3_r3] - _b[enroll_sup:`s'.sex#c.terwi2_r3])
	local b_`s' = `=round(`r(estimate)', 0.01)'
	local se_`s' = `=round(`r(se)', 0.01)'
	local pval_`s' = `r(p)'
	
	local c = 1
	foreach coef in b se { // coef then se
		mat ModelC_wealth[`=5+`s'*2+`c'',1] = ``coef'_`s''
		local c = 2
	}

	mat ModelC_wealth[`=5+`s'*2+1',2] = `=cond(`pval_`s''<0.01, 3, `=cond(`pval_`s''<0.05, 2, `=cond(`pval_`s''<0.1, 1, 0)')')'

}

qui test _b[enroll_sup:1.sex#c.terwi2_r3] = _b[enroll_sup:0.sex#c.terwi2_r3]
mat ModelC_wealth[5,1] = `=round(`r(p)', 0.001)'

qui test (_b[enroll_sup:1.sex#c.terwi3_r3] - _b[enroll_sup:1.sex#c.terwi2_r3]) = (_b[enroll_sup:0.sex#c.terwi3_r3] - _b[enroll_sup:0.sex#c.terwi2_r3])
mat ModelC_wealth[10,1] = `=round(`r(p)', 0.001)'

}

** testing effects of skills
{
qui test _b[enroll_sup:1.sex#skills_cr_r3] = _b[enroll_sup:0.sex#skills_cr_r3]
local p_skills_cr = `r(p)'
qui test _b[enroll_sup:1.sex#skills_nc2_r3] = _b[enroll_sup:0.sex#skills_nc2_r3]
local p_skills_nc2 = `r(p)'
}

** storing skill production technology
{
cap mat drop SkillProd
mat SkillProd = J(10,20,.)

forvalues s = 0/1 {
	
	* R1

	local v = 0
	foreach var in terwi2_r1 terwi3_r1 {
	
		qui lincom _b[skills_c_r1:`s'.sex#c.`var']
		local b = `r(estimate)'
		local se = `r(se)'
		local pval = `r(p)'
		
		local c = 0
		foreach coef in b se {
			mat SkillProd[`=7+`v'*2+`c'',`=1+`s'*2'] = `=round(``coef'', 0.01)'
			local ++c
		}
	
		mat SkillProd[`=7+`v'*2',`=2+`s'*2'] = `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'
	
		local ++v
	}
	
	local t = 0
	foreach type in cr nc2 {
	
		* R2
		
		local v = 0
		foreach var in skills_c_r1 parinv_r1 terwi2_r2 terwi3_r2 {
		
			if "`var'" == "parinv_r1" {
				local ++v
			}
			
			qui lincom _b[skills_`type'_r2:`s'.sex#c.`var']
			local b = `r(estimate)'
			local se = `r(se)'
			local pval = `r(p)'
			
			local c = 0
			foreach coef in b se {
				mat SkillProd[`=1+`v'*2+`c'',`=5+`t'*4+`s'*2'] = `=round(``coef'', 0.01)'
				local ++c
			}
		
			mat SkillProd[`=1+`v'*2',`=6+`t'*4+`s'*2'] = `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'
		
			local ++v
		}
		
		* R3

		local v = 0
		foreach var in skills_cr_r2 skills_nc2_r2 parinv_r2 terwi2_r3 terwi3_r3 {
		
			qui lincom _b[skills_`type'_r3:`s'.sex#c.`var']
			local b = `r(estimate)'
			local se = `r(se)'
			local pval = `r(p)'
			
			local c = 0
			foreach coef in b se {
				mat SkillProd[`=1+`v'*2+`c'',`=13+`t'*4+`s'*2'] = `=round(``coef'', 0.01)'
				local ++c
			}
		
			mat SkillProd[`=1+`v'*2',`=14+`t'*4+`s'*2'] = `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'
		
			local ++v
		}
		
		local ++t		
	}	
}	
	
}

*** Outputting results

* Main results
{
qui putexcel set "$output/Table-MainResults.xlsx", replace

putexcel C4 = matrix(ModelA_main), hcenter
putexcel G4 = matrix(ModelB_main), hcenter
putexcel K4 = matrix(ModelC_main), hcenter

qui putexcel 	C2 = "A" 		///
				G2 = "B" 		///
				K2 = "C" 		///	
				C3 = "Male"		///
				E3 = "Female"	///
				G3 = "Male"		///
				I3 = "Female"	///
				K3 = "Male"		///
				M3 = "Female"	///
				, hcenter

local v = 0
foreach var in terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 migrated_r1 indigena_r1 order_1 msib_r3 fsib_r3 skills_cr_r3 skills_nc2_r3 {
	qui putexcel B`=4+`v'*2' = "`var'"
	local ++v	
}
}

* Wealth analysis
{
qui putexcel set "$output/Table-WealthAnalysis.xlsx", replace

qui putexcel C3 = matrix(ModelA_wealth), hcenter
qui putexcel E3 = matrix(ModelB_wealth), hcenter
qui putexcel G3 = matrix(ModelC_wealth), hcenter

qui putexcel 	C2 = "A"	///
				E2 = "B"	///
				G2 = "C"	///
				, hcenter

qui putexcel 	B3 = "ter2 - m"			///
				B5 = "ter2 - f"			///
				B7 = "pval m = f"		///
				B8 = "ter3 - ter2 - m"	///
				B10 = "ter3 - ter2 - f"	///
				B12 = "pval m = f"		
}
	
* Skill production
{
qui putexcel set "$output/Table-SkillProd.xlsx", replace

putexcel C4 = matrix(SkillProd), hcenter

qui putexcel 	C2 = "Cognitive skills - R1"		///
				G2 = "Cognitive skills - R2"		///
				K2 = "Non-cognitive skills - R2"	///
				O2 = "Cognitive skills - R3"		///
				S2 = "Non-cognitive skills - R3"	
				
qui putexcel	C3 = "male"							///
				E3 = "female"						///
				G3 = "male"							///
				I3 = "female"						///
				K3 = "male"							///
				M3 = "female"						///
				O3 = "male"							///
				Q3 = "female"						///
				S3 = "male"							///
				U3 = "female"						///
				, hcenter

qui putexcel 	B4 = "cskills"	///
				B6 = "ncskills"	///
				B8 = "parinv"	///
				B10 = "terwi2"	///
				B12 = "terwi3"	
}

*===============================================================================
* Mechanisms
*===============================================================================

*** credit contraints

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(creditcons2 terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop MechCredit
mat MechCredit = J(18,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 creditcons2 {
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat MechCredit[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat MechCredit[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

*** Family disadvantage

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(parents teenagemom partnerdrunk terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop MechFamily
mat MechFamily = J(18,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 parents teenagemom partnerdrunk {

	if "`var'" == "parents" {
		local ++i
	}
		
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat MechFamily[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat MechFamily[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

*** distance (continuous)

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 zhfa_r1 indigena_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 zhfa_r1 indigena_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(educdist_r3 terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop MechDist
mat MechDist = J(18,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 educdist_r3 {
	
	if "`var'" == "educdist_r3" {
	    local i = `i' + 4
	}
		
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat MechDist[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat MechDist[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

*** All

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(parents teenagemom partnerdrunk creditcons2 educdist_r3 terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop MechAll
mat MechAll = J(18,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 creditcons2 parents teenagemom partnerdrunk educdist_r3 {
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat MechAll[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat MechAll[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

*** Outputting results

{
qui putexcel set "$output/Table-Mechanisms.xlsx", replace

putexcel C4 = matrix(MechCredit), hcenter
putexcel G4 = matrix(MechFamily), hcenter
putexcel K4 = matrix(MechDist), hcenter
putexcel O4 = matrix(MechAll), hcenter

qui putexcel 	C2 = "A" 		///
				G2 = "B" 		///
				K2 = "C" 		///	
				O2 = "D" 		///	
				C3 = "Male"		///
				E3 = "Female"	///
				G3 = "Male"		///
				I3 = "Female"	///
				K3 = "Male"		///
				M3 = "Female"	///
				O3 = "Male" 	///
				Q3 = "Female"	///
				, hcenter

local v = 0
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 creditcons2 parents teenagemom partnerdrunk educdist_r3 {
	qui putexcel B`=4+`v'*2' = "`var'"
	local ++v	
}
}

*===============================================================================